choroid Plexus

In this part of the analysis we apply Revelio algorithm to explore cell cycle dynamic of pallial and hem domain radial glial cells

Load libraries

library(Seurat)
library(Revelio)
library(princurve)
library(orthologsBioMART)
library(Matrix)
library(dplyr)
library(RColorBrewer)
library(ggplot2)
library(ggExtra)
library(cowplot)
library(wesanderson)

#Set ggplot theme as classic
theme_set(theme_classic())

Load and filter progenitors data

Hem.data <- readRDS("../QC.filtered.cells.RDS")

DimPlot(object = Hem.data,
        group.by = "Cell_ident",
        reduction = "spring",
        cols = c("#83c3b8", #"ChP"
                 "#009fda", #"ChP_progenitors"
                 "#68b041", #"Dorso-Medial_pallium"
                 "#e46b6b", #"Hem"
                 "#e3c148", #"Medial_pallium"
                 "#b7d174", #2
                 "grey40", #4
                 "black", #5
                 "#3e69ac" #"Thalamic_eminence"
                 ))

Idents(Hem.data) <- Hem.data$Cell_ident

Fit Pseudotime axis on ChP cells

ChP.data <-  subset(Hem.data, idents = c("ChP", "ChP_progenitors"))

DimPlot(ChP.data,
        reduction = "spring",
        pt.size = 1,
        cols =  c("#83c3b8", "#009fda")) + NoAxes()

Exclude septal cells

FeaturePlot(object = ChP.data ,
            features = c("Fgf8", "Fgf17", "Adamts15", "Fgfbp1"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes() & NoLegend()

ChP.data <- AddModuleScore(ChP.data,
                           features = list(c("Fgf8", "Fgf17", "Adamts15", "Fgfbp1")),
                           ctrl = 10,
                           name = "Septum")

FeaturePlot(object = ChP.data ,
            features = c("Septum1"),
            pt.size = 0.5,
            cols = rev(brewer.pal(10,"Spectral")),
            reduction = "spring",
            order = T) & NoAxes()

hist(ChP.data$Septum1, breaks = 20)

ChP.data$Septal.prog <- ChP.data$Septum1 > 0.1

p1 <- DimPlot(ChP.data,
        reduction = "spring",
        group.by = "Septal.prog",
        pt.size = 1) + NoAxes()

p2 <- FeaturePlot(object = ChP.data ,
            features = c("Fgf17"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes() & NoLegend()

p1 + p2

ChP.data <- subset(ChP.data,
                   subset = Septal.prog == FALSE & ChP.data$Spring_1 > 1300)

DimPlot(ChP.data,
        reduction = "spring",
        pt.size = 1,
        cols =  c("#83c3b8", "#009fda")) + NoAxes()

## Fit principal curve

Trajectories.ChP <- ChP.data@meta.data %>%
                    select("Barcodes", "nUMI", "Spring_1", "Spring_2")

fit <- principal_curve(as.matrix(Trajectories.ChP[,c("Spring_1", "Spring_2")]),
                       smoother='lowess',
                       trace=TRUE,
                       f = 0.8, 
                       stretch=2)

## Starting curve---distance^2: 137968435096
## Iteration 1---distance^2: 49449688
## Iteration 2---distance^2: 45168261
## Iteration 3---distance^2: 44040505
## Iteration 4---distance^2: 43717320
## Iteration 5---distance^2: 43609264
## Iteration 6---distance^2: 43575985

#The principal curve smoothed
ChP.pc.line <- as.data.frame(fit$s[order(fit$lambda),]) 

#Pseudotime score
Trajectories.ChP$Pseudotime <- fit$lambda/max(fit$lambda)

#Inverse the score if positive correlation with progenitor marker
if (cor(Trajectories.ChP$Pseudotime, ChP.data@assays$SCT@data['Hmga2', Trajectories.ChP$Barcodes]) > 0) {
  Trajectories.ChP$Pseudotime <- -(Trajectories.ChP$Pseudotime - max(Trajectories.ChP$Pseudotime))
}

ChP.data$Pseudotime <- Trajectories.ChP$Pseudotime

FeaturePlot(object = ChP.data,
            features = "Pseudotime",
            pt.size = 2,
            cols = rev(colorRampPalette(brewer.pal(n =10, name = "Spectral"))(100)),
            reduction = "spring",
            order = T) & NoAxes()

Subset progenitors and fit cell cycle axis

Prog.data <-  subset(ChP.data, idents = c("ChP_progenitors"))

DimPlot(Prog.data,
        reduction = "spring",
        pt.size = 1,
        cols =  c("#009fda")) + NoAxes()

Prog.data <- NormalizeData(Prog.data, normalization.method = "LogNormalize", scale.factor = 10000, assay = "RNA")

Prepare data for revelio input

Find mous ortologues to provided human cell cycle genes

Cellcyclegenes <- revelioTestData_cyclicGenes
head(Cellcyclegenes)

## # A tibble: 6 × 5
##   G1.S    S        G2        G2.M    M.G1  
##   <fct>   <fct>    <fct>     <fct>   <fct> 
## 1 ABCA7   ABCC2    ALKBH1    ADH4    AFAP1 
## 2 ACD     ABCC5    ANLN      AHI1    AGFG1 
## 3 ACYP1   ABHD10   AP3D1     AKIRIN2 AGPAT3
## 4 ADAMTS1 ACPP     ARHGAP11B ANKRD40 AKAP13
## 5 ADCK2   ADAM22   ARHGAP19  ANLN    AMD1  
## 6 ADCY6   ANKRD18A ARL4A     ANP32B  ANP32E

We use orthologsBioMART library to map human to mouse mouse orthologs

G1.S <- findOrthologsHsMm(from_filters = "hgnc_symbol",
                          from_values = as.character(Cellcyclegenes$G1.S), 
                          to_attributes = "external_gene_name")

S <- findOrthologsHsMm(from_filters = "hgnc_symbol",
                          from_values = as.character(Cellcyclegenes$S), 
                          to_attributes = "external_gene_name")

G2 <- findOrthologsHsMm(from_filters = "hgnc_symbol",
                          from_values = as.character(Cellcyclegenes$G2), 
                          to_attributes = "external_gene_name")

G2.M <- findOrthologsHsMm(from_filters = "hgnc_symbol",
                          from_values = as.character(Cellcyclegenes$G2.M), 
                          to_attributes = "external_gene_name")

M.G1 <- findOrthologsHsMm(from_filters = "hgnc_symbol",
                          from_values = as.character(Cellcyclegenes$M.G1), 
                          to_attributes = "external_gene_name")


gene.list <- list(G1.S$external_gene_name,
                  S$external_gene_name,
                  G2$external_gene_name,
                  G2.M$external_gene_name,
                  M.G1$external_gene_name)

CCgenes <- t(plyr::ldply(gene.list, rbind))

colnames(CCgenes) <- colnames(Cellcyclegenes)

Export counts matrix

rawCounts <- as.matrix(Prog.data[["RNA"]]@counts)

# Filter genes expressed by less than 10 cells
num.cells <- Matrix::rowSums(rawCounts > 0)
genes.use <- names(x = num.cells[which(x = num.cells >= 10)])
rawCounts <- rawCounts[genes.use, ]

rm(list = ls()[!ls() %in% c("rawCounts", "CCgenes", "ChP.data", "Prog.data")])
gc()

##             used  (Mb) gc trigger   (Mb)  max used   (Mb)
## Ncells   5781851 308.8   10748997  574.1   8259915  441.2
## Vcells 119685790 913.2  611435808 4664.9 621902469 4744.8

Run Revelio

We can now follow the tutorial form the package github page

myData <- createRevelioObject(rawData = rawCounts,
                              cyclicGenes = CCgenes,
                              lowernGeneCutoff = 0,
                              uppernUMICutoff = Inf,
                              ccPhaseAssignBasedOnIndividualBatches = F)

## 2021-12-09 17:27:00: reading data: 2.32secs

rm("rawCounts")
gc()

##             used  (Mb) gc trigger   (Mb)  max used   (Mb)
## Ncells   5783930 308.9   10748997  574.1   8259915  441.2
## Vcells 119738779 913.6  489148647 3732.0 621902469 4744.8

The getCellCyclePhaseAssignInformation filter “outlier” cells for cell cycle phase assignation. We modified the function to keep all cells as we observed this does not affect the global cell cycle fitting procedure

source("../Functions/functions_InitializationCCPhaseAssignFiltering.R")

myData <- getCellCyclePhaseAssign_allcells(myData)

## 2021-12-09 17:27:07: assigning cell cycle phases: 13.65secs

myData <- getPCAData(dataList = myData)

## 2021-12-09 17:27:29: calculating PCA: 10.25secs

myData <- getOptimalRotation(dataList = myData)

## 2021-12-09 17:27:56: calculating optimal rotation: 1.47secs

gc()

##             used   (Mb) gc trigger   (Mb)  max used   (Mb)
## Ncells   5813634  310.5   10748997  574.1  10748997  574.1
## Vcells 161090841 1229.1  489148647 3732.0 621902469 4744.8

Compare results between different ordering approched applied to the DC space

Comparison of Revelio cell cycle ordering and Angle method

To obtain the best linear ordering of cell along G1 to M phase we compare the output `ccPercentageUniformlySpaced’ form Revelio and the one obtained with the Angle method describe in Saelens et. al. Nature Biotechnology 2019

New dataframe

CellCycledata <- cbind(as.data.frame(t(myData@transformedData$dc$data[1:2,])),
                       nUMI= myData@cellInfo$nUMI,
                       Revelio.phase = factor(myData@cellInfo$ccPhase, levels = c("G1.S", "S", "G2", "G2.M", "M.G1")),
                       Revelio.cc= myData@cellInfo$ccPercentageUniformlySpaced,
                       Seurat.cc= Prog.data@meta.data[myData@cellInfo$cellID,"CC.Difference"])

Compute cycling trajectory with Angle

CellCycledata$Angle.cc <- atan2(CellCycledata$DC1, CellCycledata$DC2) / 2 / pi + .5

Cells distribution in the DC1-DC2 space

ggplot(CellCycledata, aes(DC1, DC2)) +
        geom_point(aes(color = Revelio.phase)) +
        scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5]))

p1 <- ggplot(CellCycledata, aes(DC1, DC2)) +
        geom_point(aes(color = Revelio.phase)) +
        scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5]))

p2 <- ggplot(CellCycledata, aes(DC1, DC2)) +
        geom_point(aes(color=Seurat.cc), size=2, shape=16) + 
        scale_color_gradientn(colours=rev(colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)),
                              name='Seurat_cc')

p3 <- ggplot(CellCycledata, aes(DC1, DC2)) +
        geom_point(aes(color=Angle.cc), size=2, shape=16) + 
        scale_color_gradientn(colours=rev(colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)),
                              name='Angle.cc')

p4 <- ggplot(CellCycledata, aes(DC1, DC2)) +
        geom_point(aes(color=Revelio.cc), size=2, shape=16) + 
        scale_color_gradientn(colours=rev(colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)),
                              name='Revelio_cc')


(p1+p2)/(p3+p4)

ggplot(CellCycledata, aes(x= Revelio.cc, y= nUMI/10000)) +
        geom_point(aes(color= Revelio.phase), size=0.5) +
        scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) +
        geom_smooth(method="loess", n= 50, fill="grey") +
        ylim(0,NA)

Import coordinates

Prog.data$Revelio.DC1 <- CellCycledata$DC1
Prog.data$Revelio.DC2 <- CellCycledata$DC2

Prog.data$Revelio.phase <- CellCycledata$Revelio.phase

Prog.data$Angle.cc <- CellCycledata$Angle.cc
Prog.data$Revelio.cc <- CellCycledata$Revelio.cc

p1 <- FeaturePlot(object = Prog.data,
            features = "Revelio.cc",
            pt.size = 1,
            cols = rev(brewer.pal(10,"Spectral")),
            reduction = "spring",
            order = T) & NoAxes()

p2 <- DimPlot(object = Prog.data,
        group.by = "Revelio.phase",
        pt.size = 1,
        reduction = "spring",
        cols =  c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) & NoAxes()

p3 <- FeaturePlot(object = Prog.data,
            features = "Pseudotime",
            pt.size = 2,
            cols = rev(colorRampPalette(brewer.pal(n =10, name = "Spectral"))(100)),
            reduction = "spring",
            order = T) & NoAxes()

p1 + p2 + p3 

Trajectories.progenitors <- Prog.data@meta.data %>%
                              select(Barcodes, nUMI, Spring_1, Spring_2, Pseudotime) %>% 
                              mutate(Cycling.axis= Prog.data$Revelio.cc,
                                     Phase = Prog.data$Revelio.phase,
                                     Gmnc= Prog.data@assays$RNA@data["Gmnc",],
                                     Ttr= Prog.data@assays$RNA@data["Ttr",],
                                     Htr2c= Prog.data@assays$RNA@data["Htr2c",],
                                     Top2a= Prog.data@assays$RNA@data["Top2a",])

p1 <- ggplot(Trajectories.progenitors, aes(x= Pseudotime, y= Cycling.axis)) +
        geom_point(aes(color= Phase), size=1.5) +
        scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5]))

p2 <- Trajectories.progenitors %>% arrange(Gmnc) %>%
      ggplot(aes(x= Pseudotime, y= Cycling.axis)) +
        geom_point(aes(color=Gmnc), size=1.5) +
        scale_color_gradientn(colours =c("grey90", brewer.pal(9,"YlGnBu")))

p3 <- Trajectories.progenitors %>% arrange(Ttr) %>%
      ggplot(aes(x= Pseudotime, y= Cycling.axis)) +
        geom_point(aes(color=Ttr), size=1.5) +
        scale_color_gradientn(colours =c("grey90", brewer.pal(9,"YlGnBu")))

p4 <- Trajectories.progenitors %>% arrange(Htr2c) %>%
      ggplot(aes(x= Pseudotime, y= Cycling.axis)) +
        geom_point(aes(color=Htr2c), size=1.5) +
        scale_color_gradientn(colours =c("grey90", brewer.pal(9,"YlGnBu")))

p1 + p2 + p3 + p4  + patchwork::plot_layout(ncol = 2)

p1 <- ggplot(Trajectories.progenitors, aes(x= Cycling.axis, y= Gmnc)) +
        geom_point(aes(color= Phase), size=0.5) +
        scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) +
        geom_smooth(method="loess", n= 50) +
        ylim(0,NA)

p2 <- ggplot(Trajectories.progenitors, aes(x= Cycling.axis, y= Top2a)) +
        geom_point(aes(color= Phase), size=0.5) +
        scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) +
        geom_smooth(method="loess", n= 50) +
        ylim(0,NA)

p1 + p2

rm(list = ls()[!ls() %in% c("Trajectories.progenitors", "ChP.data")])
gc()

##            used  (Mb) gc trigger   (Mb)  max used   (Mb)
## Ncells  5833769 311.6   10748997  574.1  10748997  574.1
## Vcells 59947732 457.4  391318918 2985.6 621902469 4744.8

Import progenitors cycling coordinates in the full dataset

ChP.data$Cycling.axis <- sapply(ChP.data$Barcodes,
                              FUN = function(x) {
                                if (x %in% Trajectories.progenitors$Barcodes) {
                                  x = Trajectories.progenitors[x, "Cycling.axis"]
                                } else {
                                  x = NA
                                  }
                              })

Session Info

#date
format(Sys.time(), "%d %B, %Y, %H,%M")

## [1] "09 December, 2021, 17,28"

#Packages used
sessionInfo()

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
## 
## Matrix products: default
## BLAS/LAPACK: /shared/ifbstor1/software/miniconda/envs/r-4.1.0/lib/libopenblasp-r0.3.15.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] wesanderson_0.3.6      cowplot_1.1.1          ggExtra_0.9           
##  [4] ggplot2_3.3.5          RColorBrewer_1.1-2     dplyr_1.0.7           
##  [7] Matrix_1.3-4           orthologsBioMART_0.1.0 data.table_1.14.2     
## [10] biomaRt_2.48.3         princurve_2.1.6        Revelio_0.1.0         
## [13] SeuratObject_4.0.3     Seurat_4.0.5          
## 
## loaded via a namespace (and not attached):
##   [1] BiocFileCache_2.0.0    plyr_1.8.6             igraph_1.2.9          
##   [4] lazyeval_0.2.2         splines_4.1.0          listenv_0.8.0         
##   [7] scattermore_0.7        GenomeInfoDb_1.30.0    digest_0.6.29         
##  [10] htmltools_0.5.2        fansi_0.5.0            magrittr_2.0.1        
##  [13] memoise_2.0.1          tensor_1.5             cluster_2.1.2         
##  [16] ROCR_1.0-11            globals_0.14.0         Biostrings_2.62.0     
##  [19] matrixStats_0.61.0     spatstat.sparse_2.0-0  prettyunits_1.1.1     
##  [22] colorspace_2.0-2       rappdirs_0.3.3         blob_1.2.2            
##  [25] ggrepel_0.9.1          xfun_0.28              crayon_1.4.2          
##  [28] RCurl_1.98-1.5         jsonlite_1.7.2         spatstat.data_2.1-0   
##  [31] survival_3.2-11        zoo_1.8-9              glue_1.5.1            
##  [34] polyclip_1.10-0        gtable_0.3.0           zlibbioc_1.40.0       
##  [37] XVector_0.34.0         leiden_0.3.9           future.apply_1.8.1    
##  [40] BiocGenerics_0.40.0    abind_1.4-5            scales_1.1.1          
##  [43] DBI_1.1.1              miniUI_0.1.1.1         Rcpp_1.0.7            
##  [46] viridisLite_0.4.0      xtable_1.8-4           progress_1.2.2        
##  [49] reticulate_1.22        spatstat.core_2.3-1    bit_4.0.4             
##  [52] stats4_4.1.0           htmlwidgets_1.5.4      httr_1.4.2            
##  [55] ellipsis_0.3.2         ica_1.0-2              farver_2.1.0          
##  [58] pkgconfig_2.0.3        XML_3.99-0.8           dbplyr_2.1.1          
##  [61] sass_0.4.0             uwot_0.1.10            deldir_1.0-6          
##  [64] utf8_1.2.2             labeling_0.4.2         tidyselect_1.1.1      
##  [67] rlang_0.4.12           reshape2_1.4.4         later_1.3.0           
##  [70] AnnotationDbi_1.54.1   munsell_0.5.0          tools_4.1.0           
##  [73] cachem_1.0.6           cli_3.1.0              generics_0.1.1        
##  [76] RSQLite_2.2.8          ggridges_0.5.3         evaluate_0.14         
##  [79] stringr_1.4.0          fastmap_1.1.0          yaml_2.2.1            
##  [82] goftest_1.2-3          knitr_1.36             bit64_4.0.5           
##  [85] fitdistrplus_1.1-6     purrr_0.3.4            RANN_2.6.1            
##  [88] KEGGREST_1.32.0        pbapply_1.5-0          future_1.23.0         
##  [91] nlme_3.1-152           mime_0.12              xml2_1.3.3            
##  [94] rstudioapi_0.13        compiler_4.1.0         filelock_1.0.2        
##  [97] curl_4.3.2             plotly_4.10.0          png_0.1-7             
## [100] spatstat.utils_2.2-0   tibble_3.1.6           bslib_0.3.1           
## [103] stringi_1.7.6          highr_0.9              lattice_0.20-44       
## [106] vctrs_0.3.8            pillar_1.6.4           lifecycle_1.0.1       
## [109] spatstat.geom_2.3-0    lmtest_0.9-39          jquerylib_0.1.4       
## [112] RcppAnnoy_0.0.19       bitops_1.0-7           irlba_2.3.3           
## [115] httpuv_1.6.3           patchwork_1.1.1        R6_2.5.1              
## [118] promises_1.2.0.1       KernSmooth_2.23-20     gridExtra_2.3         
## [121] IRanges_2.28.0         parallelly_1.29.0      codetools_0.2-18      
## [124] MASS_7.3-54            assertthat_0.2.1       withr_2.4.3           
## [127] sctransform_0.3.2      S4Vectors_0.32.0       GenomeInfoDbData_1.2.6
## [130] mgcv_1.8-36            parallel_4.1.0         hms_1.1.1             
## [133] grid_4.1.0             rpart_4.1-15           tidyr_1.1.4           
## [136] rmarkdown_2.11         Rtsne_0.15             Biobase_2.52.0        
## [139] shiny_1.7.1

---

1. Institute of Psychiatry and Neuroscience of Paris, INSERM U1266, 75014, Paris, France, matthieu.moreau@inserm.fr↩︎